Method and apparatus for ultra fast symmetry and SIMD based projection-back projection for 3D pet image reconstruction

ABSTRACT

A method and apparatus are provided for reconstructing 3D image. The method may include the steps of: detecting a plurality of line of response (LORs) emitted from an object; transforming the plurality of LORs into first sinogram data; back-projecting the first sinogram data with a plurality of projection angles to produce image data for the object; and projecting the produced image data with the plurality of projection angles to transform the image data into second sinogram data. The back-projecting may include filling pixels of image plane for each of the plurality of projection angles with the first sinogram data and rotating a coordinate axis of the image plane with corresponding projection angle to produce the image data. The projecting may include rotating the image data with corresponding projection angle in an opposite direction before projecting the image data with the plurality of projection angles. The projecting and the back-projecting may use symmetry properties in coordinate space.

The present application claims priority from Korean Patent ApplicationNo. 10-2006-0042155 filed on May 10, 2006, Korean Patent Application No.10-2007-0027305 filed on Mar. 20, 2007, and Korean Patent ApplicationNo. 10-2007-0040515 filed on Apr. 25, 2007, the entire subject matter ofwhich is incorporated herein by reference.

BACKGROUND

1. Field

The present invention relates to a method and apparatus for ultra fastsymmetry and SIMD based projection-backprojection for 3D PET imagereconstruction. More specifically, the present invention relates to amethod and apparatus based on the symmetry properties of the projectionand backprojection processes, especially in the 3D OSEM algorithm thatrequires a plurality of projections and back-projections.

2. Background

There has been a remarkable progress in PET development over the recentyears, especially in the areas of hardware, software and computerimplementation of image reconstruction. Recent developments in PETscanners (e.g., HRRT (High Resolution Research Tomograph) developed byCTI (now Siemens)) allow greatly enhanced resolution and sensitivity. Insuch PET scanners, the amount of collected coincidence line datacontains more than 4.5×10⁹ coincidence lines of response generated by asmany as 120,000 nuclear detectors. Such large amount of data and thereconstruction of this data set pose to be a real problem in HRRT. Thatis, they pose to be major problems in achieving further developments andapplications of high resolution PET scanners. Thus, in such types ofscanners, obtaining one set of reconstructed images often requires manyhours of image reconstruction. For example, in HRRT with full datacollection in normal brain scans (using SPAN 3), the imagereconstruction time is almost eighty minutes. This makes it practicallyimpossible to attempt any list mode based dynamic imaging since theimage reconstruction time takes many days (as long as 43 hours or morefor 32 frame dynamic image reconstruction).

In general, tomographic images can be reconstructed by two approaches,one being an analytic method and the other being an iterative approach.PET scanners of different types were developed in the mid 1970's and theapplication of various tomographic image reconstruction techniques wasnaturally introduced in the field. In case of an analytic approach suchas backprojection and filtering or filtered back-projection (FB), anartifact known as a streak artifact is frequently generated. This isespecially true when the detector arrangements are not uniform such asin the case of HRRT (e.g., Siemens' High Resolution Research Tomograph)where the detectors are arranged in a set of blocks in an octagonalshape. These types of detector arrangements often involve missing datadue to the gaps between the blocks and result in a severe streakartifact in the case of the FB technique. Therefore, alternativeapproaches such as an EM (Expectation Maximum) algorithm have beensought. Generally, the EM approach requires several steps in thereconstruction process, the two major steps of which are: projection(forward projection) to create projection data from the image or objectand backprojection into the image domain for the final imagereconstruction. In the EM algorithms, these two processes are repeateduntil satisfactory images are obtained. Obviously, these repeatedprojection and backprojection processes are time consuming and have beenthe major drawback of the EM approach compared to straight filteredbackprojection (FB) algorithm. In addition, in case of 3D imagereconstruction, the computational burden increases out of proportion dueto the astronomical increases in the coincidence lines or the line ofresponses (LOR). This is a major stumbling block in the daily operationof high resolution PET scanners. Thus, there is a strong need forimproving the computational speed or the reconstruction time in EMapproaches, especially with high end PET scanners such as HRRT.

Projection methods usually employ a system matrix, which is determinedby the geometric factor of the scanner. As the resolution of the PETimage improves and the number of slices increases, the size of thematrix is also increased drastically in proportion to the increases inLORs, thus resulting in not only the need for a large memory but alsothe total computation time. Current HRRT, for example, requires nearlyeighty minutes of reconstruction time, in addition to the generation ofsinograms and appropriate data streaming processes such as attenuation,random and scatter corrections, a set of precursors to thereconstruction processes. To remedy the computational burden of imagereconstruction, a number of alternative proposals such as linearintegration have been proposed, as well as the use of multiple CPUs or acluster computer system approach. Most of the techniques, however, arenot practically useful since such cluster computing requires a largedata transfer time, although the overall computation is faster than asingle unit.

Obtaining or generating projection data can be divided into twocategories, namely, ray-driven method and pixel-driven method. Theray-driven method calculates the linear integration directly along theray path connecting the centers of the two opposite detector cells,whereas the pixel-driven method calculates the linear integration alongthe ray path centered around an image pixel for the entire projectionangles. The ray-driven method is often used in projection, whilepixel-driven method is used in backprojection.

In early reconstruction techniques, projection was obtained by weighingthe ray passing through the areas of pixels with the assumption that theray path is a strip with a finite width. It, however, involves a largeamount of computation as well as the storage of a large number of matrixor data. Concurrently, Shepp and Logan proposed a simple andcomputationally efficient algorithm, which requires computing the lengthof the ray path intersecting with each pixel (instead of the areas).

To speed up the computation, there have been a number of attempts toreduce the reconstruction time. An incremental algorithm has also beendeveloped in which the symmetric property between the neighboring pixelsis considered to calculate the position of intersection of a ray. Thisidea was expanded to 3D reconstruction in cylindrical geometry usingoblique rays. In 3D form with a multi-ring system such as HRRT, itbecame apparent that true 3D approaches will be required to fullyutilize the oblique rays to thereby improve the statistics of the image.

BRIEF DESCRIPTION OF THE DRAWINGS

Arrangements and embodiments may be described in detail with referenceto the following drawings in which like reference numerals refer to likeelements and wherein:

FIG. 1A illustrates a 3-D object and its projection to a 2-D projectionplane.

FIG. 1B illustrates a y′-z plane view of FIG. 1A.

FIG. 1C illustrates an example of the transverse plane at θ=0° and theline integrals along the y′ line projected on to the x_(r).

FIG. 2A illustrates a rotation of projection ray (ray-path) or frameonto the image plane which is on the fixed reference (x,y) coordinate.

FIG. 2B illustrates the case where the projection ray (ray-path) framecoincides with the fixed (x′,y′) coordinate.

FIG. 3 illustrates the relationships among z, y′, x′, θ, x_(r), y_(r),r_(y) _(r) _(,n,θ), Z_(y) _(r) _(,n,θ) by using the y′-z plane view.

FIG. 4A illustrates an example of the “Mirror-symmetry” used in theproposed SSP method.

FIG. 4B illustrates an example of the “φ-symmetry” used in the proposedSSP method.

FIG. 5A illustrates the y′-symmetry in the proposed SSP method.

FIG. 5B illustrates the θ-symmetry in the proposed SSP method.

FIG. 6 illustrates the concept of the balanced job distribution based onthe sum of the ray path length.

FIG. 7 is a flow chart of the projection according to an embodiment ofthe present invention.

FIG. 8 is a flow chart of the back-projection according to an embodimentof the present invention.

FIG. 9A illustrates a comparison of projection data between the existingmethod and the proposed SSP method.

FIG. 9B illustrates cut-views of sinogram at a specific view.

FIG. 10A illustrates a comparison of simple back-projection imagesbetween the existing method and the proposed SSP method.

FIG. 10B illustrates cut-views (profiles) at an x axis (y=154, z=103).

FIG. 11A illustrates a set of reconstructed images with the existingmethod, the proposed SSP method and the differences.

FIG. 11B illustrates cut-views (profiles) at an x axis (y=154, z=103).

FIG. 12 illustrates the symmetry relationship in SIMD.

DETAILED DESCRIPTION

A detailed description may be provided with reference to theaccompanying drawings. One of ordinary skill in the art may realize thatthe following description is illustrative only and is not in any waylimiting. Other embodiments of the present invention may readily suggestthemselves to such skilled persons having the benefit of thisdisclosure.

Overview of the Projection, Backprojection and Symmetry Properties

A. Overview of Projection and Backprojection in the Aligned (Reference)Frame with Rotated Projection Plane

FIG. 1A illustrates a 3-D object and its projection to a 2-D projectionplane. FIG. 1B is a y′-z plane view of FIG. 1A. There is provided arelation between projection planes and the path of the projection ray,i.e., {right arrow over (ι)}_(x) _(r) _(,y) _(r) _(,φ,θ). The lineintegral will be performed along {right arrow over (ι)}_(x) _(r) _(,y)_(r) _(,φ,θ). The angle +θ indicates the oblique angle of the imageplanes toward an upper side. x_(r) and y_(r) are the coordinates of the2-D projection plane of the 3-D object. FIG. 1C is an example of thetransverse plane at θ=0° and the line integrals along the y′ lineprojected on to the x_(r). In this figure, φ indicates rotation of theprojection ray set in reference to the coordinate axis (x,y). Dottedlines with arrow indicate the projection rays. Each projection ray isdetermined by four variables, i.e., (x_(r),y_(r),φ,θ). The projectionplane shown in FIG. 1A is composed of 2D projection data (or bed ofnails) in coordinate (x_(r),y_(r)) at a given θ and φ.

In 3D tomographic image processing, the projection ray has 4 dimensions,namely, x_(r), y_(r), φ and θ, as shown in the FIG. 1. Projection is aprocess of converting a 3D object function or image information in 3Dcoordinates to 2D projection coordinates, and the projection plane asshown in FIG. 1.

With a fixed view angle φ and oblique angle θ, the projection plane isdetermined as shown in FIG. 1A. In this projection plane, the definitionof the projection operation can be expressed in the form of linearintegral given by: $\begin{matrix}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {\int{\int{\int{{I\left( {x,y,z} \right)}{\delta\left( {\overset{\longrightarrow}{\iota}}_{x_{r},y_{r},\phi,\theta} \right)}{\mathbb{d}x}{\mathbb{d}y}{\mathbb{d}z}}}}}} & (1)\end{matrix}$where

{right arrow over (ι)}_(x) _(r) _(,y) _(r) _(,φ,θ) is the ray path,

x_(r)=x′=x cos φ+y sin φ,

y_(r)=z−(−x sin φ+y cos φ)tan θ.

Equation (1) describes that projection P_(φ,θ)(x_(r),y_(r)) is a sum ofthe pixel values in an image function I(x,y,z) along the path ofprojection ray, {right arrow over (ι)}_(x) _(r) _(,y) _(r) _(,φ,θ), inthe image domain. Here, δ(.) represents a sampling function.

FIG. 2 illustrates the concepts of the rotating projection ray(ray-path) frame onto the image plane in (x,y) and the proposed fixed(aligned) projection frame with a rotating image plane. FIG. 2Aillustrates the rotation of projection ray (ray-path) or frame onto theimage plane, which is on the fixed reference (x,y) coordinate. This isthe conventional scheme applied to most of the image reconstructions.FIG. 2B illustrates a case where the projection ray (ray-path) framecoincides with the fixed (x′,y′) coordinate. The latter means that theimage plane is now rotated instead of the reference frame, theprojection ray frame. This latter scheme is the basis of the proposedSSP (Symmetry and SIMD based Projection-backprojection) method. Thisscheme allows the symmetry property of the image plane to be utilized,thereby reducing the overall projection and backprojection time. In caseof backprojection, the image function I (x′,y′,z) represents theintermediate stage of an image to be reconstructed and is a temporaryimage data.

It is known that a projection plane at angle θ can be rotated (oraligned) against reference axes, (x,y). Further, x′ coincides with x_(r)in FIG. 2. The well known relation between the rotated coordinate(x′,y′,z) and the image coordinate (x,y,z) in the cylindrical coordinateare provided by the following: $\begin{matrix}{\begin{pmatrix}x^{\prime} \\y^{\prime} \\z\end{pmatrix} = {{R_{\phi} \cdot \begin{pmatrix}x \\y \\z\end{pmatrix}} = {\begin{pmatrix}{\cos\quad\phi} & {\sin\quad\phi} & 0 \\{{- \sin}\quad\phi} & {\cos\quad\phi} & 0 \\0 & 0 & 1\end{pmatrix}\begin{pmatrix}x \\y \\z\end{pmatrix}}}} & (2)\end{matrix}$where

R_(φ) is the rotation matrix.

Equation (1) can then be simplified as the weighed sum along the raypath and can be written as follows: $\begin{matrix}\begin{matrix}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {\int{\int{\int{{I\left( {x,y,z} \right)}{\delta\left( {\overset{\longrightarrow}{\iota}}_{x_{r},y_{r},\phi,\theta} \right)}{\mathbb{d}x}{\mathbb{d}y}{\mathbb{d}z}}}}}} \\{= {\int{{I\left( {x^{\prime},y^{\prime},z} \right)}{\mathbb{d}y^{\prime}}}}}\end{matrix} & (3)\end{matrix}$where

I(x′,y′,z)=R_(φ)I(x,y,z)

x′=x_(r)

In (3),y′ is assumed to be the integration variable.

As a special case, when the oblique angle θ is equal to zero, theprojection rays are in parallel with the y′ axis (see FIG. 1). Inmulti-layer or 3D PET, oblique rays (θ≠0) are collected and used fortrue 3D reconstruction since the full utilization of the oblique rayswill enhance the image quality due to the increased statistics.Extension of (3) to a 3D case can be represented by θ≠0. It should benoted that the coordinates in the transverse plane are now independentof θ and the ray projected to the transverse plane is parallel to the y′axis. By the trigonometric relationship depicted in FIG. 3, z can bewritten as follows:z=y_(r)+y′tan θ  (4)

FIG. 3 illustrates the relationships among z, y′, x′, θ, x_(r), y_(r),r_(y) _(r) _(,n,θ) and Z_(y) _(r) _(,n,θ) by using the y′-z plane view.It should be noted that these relationships are valid on all the y′-zplanes with any x′. Also, y′ is now noted with discrete value n.

By substituting (4) into (3), the following discrete form of projectiondata can be obtained: $\begin{matrix}\begin{matrix}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {\int{{I\left( {x_{r},y^{\prime},{y_{r} + {y^{\prime}\tan\quad\theta}}} \right)}{\mathbb{d}y^{\prime}}}}} \\{\approx {\sum\limits_{n = {{- N_{I_{x_{r}}}}/2}}^{N_{I_{x_{r}}}/2}\quad\left\lbrack {I\left( {x_{r},n,{y_{r} + {n\quad\tan\quad\theta}}} \right)} \right\rbrack}} \\{\approx {\sum\limits_{n = {{- N_{I_{x_{r}}}}/2}}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{{I\left( {x_{r},n,Z_{y_{r},n,\theta}} \right)}\left( {1 - r_{y_{r},n,\theta}} \right)} +} \\{{I\left( {x_{r},n,{Z_{y_{r},n,\theta} + 1}} \right)}\left( r_{y_{r},n,\theta} \right)}\end{bmatrix}}}\end{matrix} & (5)\end{matrix}$where

z_(y) _(r) _(,n,θ)=y_(r)+n tan θ=Z_(y) _(r) _(,n,θ)+r_(y) _(r) _(,n,θ)

Z_(y) _(r) _(,n,θ)εINTEGER, r_(y) _(r) _(,n,θ)εREAL, and 0≦r_(y) _(r)_(,n,θ)<1.

Z_(y) _(r) _(,n,θ) is integer value of Z_(y) _(r) _(,n,θ)

r_(y) _(r) _(,n,θ) is interpolation coefficient or remainder of Z_(y)_(r) _(,n,θ)

r_(y) _(r) _(,n,θ)+r_(y) _(r) _(,−n,θ)=1 or r_(y) _(r) _(,−n,θ)=1−r_(y)_(r) _(,n,θ)

n is discrete value of y′ N_(I_(x_(r)))

-   -   is integration length of y′ at given y′ at given y_(r), φ, and        θ.

It should be noted that since projection rays and the coordinates(x,y,z) or (x′,y′,z) are not always coincident. Interpolation is alwaysrequired. In case of SSP implementation, a 1D linear interpolation isneeded along the z axis, as shown in (5). Note that r_(y) _(r) _(,n,θ)and 1−r_(y) _(r) _(,n,θ) are the interpolation coefficients along the gaxis.

On the other hand, the backprojection process is the transpose ofprojection. Therefore, the same idea as in the projection (or projectiondata generation) process can be applied to backprojection or image datageneration. In case of an iterative reconstruction such as the OS-EMalgorithm, a plurality of projections and backprojections are required.That is, algorithm requires many repeated image reconstructions as wellas projection data generations.

Now the backprojection from projection data, P_(φ,θ)(x_(r),y_(r)), canbe performed, as shown by the following: $\begin{matrix}\begin{matrix}{{I\left( {x,y,z} \right)} = {\int_{\phi}{\int_{\theta}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)}{\mathbb{d}\theta}{\mathbb{d}\phi}}}}} \\{= {\int_{\phi}{\int_{\theta}{P_{\phi,\theta}\left( {{{x\quad\cos\quad\phi} + {y\quad\sin\quad\phi}},} \right.}}}}\end{matrix} & (6) \\{\left. {z - {\left( {{{- x}\quad\sin\quad\phi} + {y\quad\cos\quad\phi}} \right)\tan\quad\theta}} \right){\mathbb{d}\theta}{\mathbb{d}\phi}} & \quad\end{matrix}$

By using rotation, the matrix defined in (2), (6) can also be writtenas: $\begin{matrix}{{I\left( {x,y,z} \right)} = {\int_{\phi}{\int_{\theta}{{P_{\phi,\theta}\left( {x^{\prime},{z - {y^{\prime}\quad\tan\quad\theta}}} \right)}{\mathbb{d}\theta}{\mathbb{d}\phi}}}}} & (7)\end{matrix}$

In the aligned frame, x_(r)=x′ is defined.

To simplify the notation, let us define I_(φ) as the sum of theprojection planes about θ for a fixed φ as given below as anintermediate stage image before the finally reconstructed image. Thecomputation can then be further simplified, for example, by using thesymmetry property in θ, $\begin{matrix}{{I_{\phi}\left( {x,y,z} \right)} \equiv {\int_{\theta}{{P_{\phi,\theta}\left( {x^{\prime},{z - {y^{\prime}\quad\tan\quad\theta}}} \right)}{\mathbb{d}\theta}}}} & (8)\end{matrix}$

Using the rotation matrix given in (2) once again, (8) can be noted asfollows:I_(φ)(x′,y′,z)=R_(φ)I₁₀₀ (x,y,z)or I_(φ)(x,y,z)=R_(−φ)I_(φ)(x′,y′,z)   (9)

Using (9), (7) can further be noted as: $\begin{matrix}{{I\left( {x,y,z} \right)} = {\int_{\phi}{R_{- \phi}{I_{\phi}\left( {x^{\prime},y^{\prime}\quad,z} \right)}{\mathbb{d}\phi}}}} & (10)\end{matrix}$

Equation (10) is equivalent to rotating the set of partiallyreconstructed image planes at a given angle φ and then integrating thesepartially reconstructed images over all φ's to form the finallyreconstructed image I(x,y,z).

B. Symmetry Properties

As discussed above, the full utilization of symmetry properties will bethe central part of the proposed SSP method. Hereafter, the symmetricpoints represent the similar coefficients (the same value but withdifferent polarity and/or different coordinates) that will be shared incomputation without additional calculations as well as the use ofidentical instructions. That is, the symmetric points have the sameinterpolation coefficient or complementary value. The applicant foundthat there exists sixteen symmetric points, which are practicallyusable, namely “Mirror-symmetry,” “φ-symmetry,” “y′-symmetry” and“θ-symmetry.”

(a) Mirror-symmetry (+x′ and −x′, see FIG. 4A)

One of the interesting aspects of the image reconstruction is the use ofthe symmetric properties of each image point. That is, once one point iscalculated with an appropriate interpolation coefficient, one can assignthe same coefficient value in the opposite point(s) by using symmetryproperties of various types. One of the simplest forms is the use of they′-axis symmetry, as shown in FIG. 4A. In this computation, a point,(x′₀,y′₀) at +x_(r) is identical in value (x′₁,y′₁) at −x_(r) withchange in polarity only. An illustration of this “Mirror-symmetry” isshown in FIG. 4A. In this illustration, an example is given forsymmetric points due to “Mirror-symmetry”, i.e.,(x′₀,y′₀)→(x′₁,y′₁)=(−x′₀,y′₀). This symmetry property reduces thecomputation time by half.

(b) φ-symmetry (φ and φ+90° see FIG. 4B)

Similar to the case of the “Mirror-symmetry”, one can also demonstratethat there is symmetry between the points at φ and φ+90° symmetry. Inthis case, (x′₀,y′₀)→(x′₁,y′₁)=(−y′₀, x′₀), i.e., the polarity as wellas the coordinates are changed. FIG. 4B illustrates this φ-symmetry casewith a numerical example. As shown in the figure, once a point iscalculated at φ and an interpolation coefficient value for this point isobtained, for example, at (x′₀,y′₀), the same coefficient value can beused at (x′₁,y′₁) with coordinate and polarity changes as (−y′₀,x′₀)after a 90° rotation, i.e., at φ+90°. This symmetry property once againpermits the computation time to be reduced by half.

FIG. 4 illustrates the examples of the “Mirror-symmetry” and“φ-symmetry” used in the proposed SSP method. FIG. 4A illustrates anexample of “Mirror-symmetry” against the y′ axis. That is, once acoordinate point (x′₀,y′₀) is known, a symmetric point (x′₁,y′₁) acrossthe y′ axis can be obtained or assigned without additional computation.This will speed up the computation by a factor of two or computationtime by half. FIG. 4B illustrates the symmetric property of theprojection data against 0, i.e., if an interpolation coefficient of acoordinate point (x′₀,y′₀) is known at φ, then a symmetric point(x′₁,y′₀) at φ+90° can be obtained without additional computation. Inother words, the projection data calculation at φ+90° share the samecoordinate values at φ with simple coordinate changes, i.e., (x′,y′) atφ to (−y′,x′) at φ+90°. This symmetry property reduces the computationtime by half.

(c) y′-symmetry (+n and −n, see FIG. 5A)

Unlike mirror- and φ-symmetry, y′- and θ-symmetry are in the y′-z plane.It should be noted that y′-symmetry is intimately related to θ-symmetry.Therefore, it is easy to illustrate the two together, as follows. They′-symmetry property has the effect of changing the interpolationcoefficients to complementary values. In fact, y′-symmetry has the samesymmetry property as θ-symmetry, as discussed below. First, the propertyof y′-symmetry is illustrated in FIG. 5A. As shown in FIG. 5A, variableZ can now be noted as follows (see also FIG. 3):z_(y) _(r) _(,n,θ)=y_(r)+ntanθ=Z_(y) _(r) _(,n,θ)+r_(y) _(r)_(,n,θ)  (11)wheren, Z_(y) _(r) _(,n,θ and y) _(r)εINTEGERr_(y) _(r) _(,n,θ)εREAL and 0≦r_(y) _(r) _(,n,θ)<1$\begin{matrix}\begin{matrix}{{z_{y_{r},n,\theta} + z_{y_{r},{- n},\theta}} = {\left( {y_{r} + {n\quad\tan\quad\theta}} \right) + \left( {y_{r} - {n\quad\tan\quad\theta}} \right)}} \\{= {Z_{y_{r},n,\theta} + r_{y_{r},n,\theta} + Z_{y_{r},{- n},\theta} + r_{y_{r},{- n},\theta}}} \\{= {2y_{r}\varepsilon\quad I\quad N\quad T\quad E\quad G\quad E\quad R}}\end{matrix} & (12) \\\begin{matrix}{{\therefore{r_{y_{r},n,\theta} + r_{y_{r},{- n},\theta}}} = 1} \\{\quad{r_{y_{r},{- n},\theta} = {1 - r_{y_{r},n,\theta}}}}\end{matrix} & \quad\end{matrix}$where

r_(y) _(r) _(,n,θ) and 1−r_(y) _(r) _(,n,θ) represent die interpolationcoefficients.

Using this symmetry, two interpolation operations can be simultaneouslyperformed, thereby reducing the number of loops in y′(=n) by half.Therefore, (5) can be rewritten as follows: $\begin{matrix}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {{I\left( {x_{r},0,y_{r}} \right)} + {\sum\limits_{n = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{{I\left( {x_{r},n,Z_{y_{r},n,\theta}} \right)}\left( {1 - r_{y_{r},n,\theta}} \right)} +} \\{{{I\left( {x_{r},n,{Z_{y_{r},n,\theta} + 1}} \right)}r_{y_{r},n,\theta}} +} \\{{{I\left( {x_{r},{- n},Z_{y_{r},{- n},\theta}} \right)}\left( {1 - r_{y_{r},{- n},\theta}} \right)} +} \\{{I\left( {x_{r},{- n},{Z_{y_{r},{- n},\theta} + 1}} \right)}r_{y_{r},{- n},\theta}}\end{bmatrix}}}} & (13)\end{matrix}$

Equation (13) can be transformed into a factored form to reduce thenumber of instructions required to execute the computation. Furthermore,the interpolation coefficient r_(y) _(r) _(,n,θ) can be reduced tosimple functions of only n and θ and can be derived from (11). Becausey_(r) is always an integer and r_(y) _(r) _(,n,θ) satisfies thecondition 0≦r_(y) _(r) _(,n,θ)<1, r_(y) _(r) _(,n,θ) is a function ofonly n and θ. Therefore, the notation r_(y) _(r) _(,n,θ) can besimplified to r_(n,θ) (see and compare with FIGS. 3 and 5A). Further,(1−r_(y) _(r) _(,n,θ)) is equal to r_(y) _(r) _(,−n,θ). Similarly,(1−r_(y) _(r) _(,−n,θ)) is equal to r_(n,θ). Using the simplified forms,(13) can be transformed to a factored form, as shown below:$\begin{matrix}{{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {{I\left( {x_{r},0,y_{r}} \right)} + {\sum\limits_{x = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}\left\{ {{I\left( {x_{r},n,Z_{y_{r},n,\theta}} \right)} +} \right. \\{{\left. {I\left( {x_{r},{- n},{Z_{y_{r},{- n},\theta} + 1}} \right)} \right\}\left( {1 - r_{n,\theta}} \right)} +} \\{\left\{ {{I\left( {x_{r},n,{Z_{y_{r},n,\theta} + 1}} \right)} + {I\left( {x_{r},{- n},Z_{y_{r},{- n},\theta}} \right)}} \right\} r_{n,\theta}}\end{bmatrix}}}} & (14)\end{matrix}$

(d) θ-symmetry (+θ and −θ, see FIG. 5B)

Now, let us illustrate how the θ-symmetry shares the calculation ofinterpolation coefficients with the y′-symmetry discussed above. Asshown in FIG. 5B, only the changes in the interpolation coefficientschange their own complementary values when the signs of θ change.

Using (4), (5), (11) and (12), it can be shown that z and r are givenas: $\begin{matrix}\begin{matrix}\begin{matrix}{z_{y_{r},n,\theta} = {y_{r} + {n\quad\tan\quad\theta}}} \\{= {y_{r} + {\left( {- n} \right)\quad\tan\quad\left( {- \theta} \right)}}} \\{= z_{y_{r},{- n},{- \theta}}}\end{matrix} \\\begin{matrix}{z_{y_{r},n,{- \theta}} = {y_{r} + {n\quad\tan\quad\left( {- \theta} \right)}}} \\{= {y_{r} + {\left( {- n} \right)\quad\tan\quad\theta}}} \\{= z_{y_{r},{- n},\theta}}\end{matrix}\end{matrix} & (15) \\{{Therefore},} & \quad \\\begin{matrix}{r_{{- n},\theta} = {r_{n,{- \theta}} = {1 - r_{n,\theta}}}} \\{r_{{- n},{- \theta}} = r_{n,\theta}} \\{Z_{y_{r},n,\theta} = Z_{y_{r},{- n},{- \theta}}} \\{Z_{y_{r},n,{- \theta}} = Z_{y_{r},{- n},\theta}}\end{matrix} & \quad\end{matrix}$

Using the relation given in (15), (14) can be written as follows:$\begin{matrix}{{P_{\phi,{- \theta}}\left( {x_{r},y_{r}} \right)} = {{I\left( {x_{r},0,y_{r}} \right)} + \quad{\sum\limits_{n = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{\begin{Bmatrix}{{I\left( {x_{r},n,Z_{y_{r},n,{- \theta}}} \right)} +} \\{I\left( {x_{r},{- n},{Z_{y_{r},{- n},{- \theta}} + 1}} \right)}\end{Bmatrix}\left( {1 - r_{n,{- \theta}}} \right)} +} \\{\begin{Bmatrix}{{I\left( {x_{r},n,{Z_{y_{r},n,{- \theta}} + 1}} \right)} +} \\{I\left( {x_{r},{- n},Z_{y_{r},{- n},{- \theta}}} \right)}\end{Bmatrix}r_{n,{- \theta}}}\end{bmatrix}}}} & (16) \\{\quad{= {{I\left( {x_{r},0,y_{r}} \right)} + {\sum\limits_{n = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{\begin{Bmatrix}{{I\left( {x_{r},n,Z_{y_{r},{- n},\theta}} \right)} +} \\{I\left( {x_{r},{- n},{Z_{y_{r},n,\theta} + 1}} \right)}\end{Bmatrix}r_{n,\theta}} +} \\{\begin{Bmatrix}{{I\left( {x_{r},n,{Z_{y_{r},{- n},\theta} + 1}} \right)} +} \\{I\left( {x_{r},{- n},Z_{y_{r},n,\theta}} \right)}\end{Bmatrix}\left( {1 - r_{n,\theta}} \right)}\end{bmatrix}}}}} & \quad\end{matrix}$

Equation (14) is for P_(φ,θ)(x_(r),y_(r)), while (16) is for theP_(φ,−θ) (x_(r),y_(r)). For the computation of (16), the relations givenin (15), i.e., r_(n,θ), r_(n,−θ), Z_(y) _(r) _(,n,θ), and Z_(y) _(r)_(,n,−θ) can be used once again and further simultaneously calculateP_(φ,θ)(x_(r),y_(r)) and P_(φ,−θ(x) _(r), y_(r)).

As shown above, y′- and θ-symmetries share the common properties. They′-symmetry, therefore, is used to reduce the number of loop countsy′(=n), while the θ-symmetry is used to reduce the number of loop countsθ.

In summary, a total of sixteen points of symmetry for the SSP method,namely, 2 (Mirror-symmetry), 2 (φ+90° symmetry), 2 (y′-symmetry) and 2(θ-symmetry), are used. Thus, multiplication of these leads to 16. Inaddition to these symmetry based reductions in computation time (16times), as will be detailed, SIMD will speed-up the computation furtherby another factor of 4. As such, a total speed gain of 64 (16×4) hasbeen achieved with the new SSP method.

FIG. 5 illustrates y′- and/or θ-symmetry in the proposed SSP method.FIG. 5A shows y′-symmetry, while FIG. 5B shows θ-symmetry. The left sideof each figure shows the z-y′ plane cut-view, while the right part showsan expanded view of the left. Each quarter illustrates the details ofthe interpolation coefficient, r′s. This y′- or θ-symmetry propertyindicates that if one of the interpolation coefficients is calculatedfor one point, then the remaining interpolation coefficients for theother three points can be obtained without additional computation. Theamount of calculations, therefore, can be reduced to one fourth (¼) byusing y′- and θ-symmetries.

Parallelization of Computation using SIMD for SSP method

A. SIMD (Single Instruction Multiple Data) and Symmetry Properties

It is important to reduce the number of instructions per loop forpractical purposes, while reducing the number of loops within theprojector/backprojector, to realize a fast algorithm. The number ofinstructions per loop was reduced by using the SIMD technique. The useof SIMD reduces the number of instructions per loop, thereby reducingthe overall computation time by a factor of four. The operand of SIMDshould be determined carefully to improve the efficiency. As shownabove, projection/backprojection of sixteen points were simultaneouslyperformed. For this operation, sixteen symmetry points were divided bytwo groups, one for +θ and the other for −θ. Since a single SIMDinstruction data set consists of four data sets, each group having thesame (+θ or −θ), was coupled to one of the two SIMD instructiondata-sets.

These two groups having the same θ is based on (14) and (16). FIG. 12shows these sixteen pairs of symmetry relationship between rotated imagedata sets and projection data sets. The upper group (a-1 and a-2) sharesthe same θ, while the lower group shares −θ. Each group is divided bytwo sub-groups. Each upper subgroup shares the same linear-interpolationcoefficients, while each lower group shares the complementary values tothat of the each upper group. Each single subgroup, therefore, can behandled as a single SIMD data set (or single SLMD pack). Z_(y) _(r)_(,y′,θ) and z_(y) _(r) _(,−y′,−θ) share the same value. Similarly,z_(y) _(r) _(,y′,θ) and z_(y) _(r) _(,−y′,θ) also share the same value.

As noted above, sixteen symmetry pairs are divided into two groups, eachsharing the same θ. In addition, these two groups were once againdivided into two subgroups. The members of subgroups a-1 and b-2 sharethe same interpolation coefficient (e.g., r_(y′,θ)), while the subgroupa-2 and b-1 have complementary values (e.g., 1−r_(y′,θ)) asinterpolation coefficients. Each subgroup represents one packed SIMDdata, as shown below: $\begin{matrix}{{{SI}_{x^{\prime}}\left( {n,z} \right)} = \left\{ {{I\left( {x^{\prime},n,z} \right)},{I\left( {{- x^{\prime}},n,z} \right)},{I\left( {{- n^{\prime}},x^{\prime},z} \right)},{I\left( {{- n},{- x^{\prime}},z} \right)}} \right\}} & (17) \\{{{SI}_{x^{\prime}}\left( {{- n},z} \right)} = \left\{ {{I\left( {x^{\prime},{- n},z} \right)},{I\left( {{- x^{\prime}},{- n},z} \right)},{I\left( {n,x^{\prime},z} \right)},{I\left( {n,{- x^{\prime}},z} \right)}} \right\}} & \quad \\{{\left( {SP} \right)_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = \left\{ {{P_{\phi,\theta}\left( {x_{r},y_{r}} \right)},{P_{\phi,\theta}\left( {{- x_{r}},y_{r}} \right)},{P_{{\phi + \frac{\pi}{2}},\theta}\left( {x_{r},y_{r}} \right)},{P_{{\phi + \frac{\pi}{2}},\theta}\left( {{- x_{r}},y_{r}} \right)},} \right\}} & \quad \\{{\left( {SP} \right)_{\phi,{- \theta}}\left( {x_{r},y_{r}} \right)} = \left\{ {{P_{\phi,{- \theta}}\left( {x_{r},y_{r}} \right)},{P_{\phi,{- \theta}}\left( {{- x_{r}},y_{r}} \right)},{P_{{\phi + \frac{\pi}{2}},\theta}\left( {x_{r},y_{r}} \right)},{P_{{\phi + \frac{\pi}{2}},\theta}\left( {{- x_{r}},y_{r}} \right)},} \right\}} & \quad\end{matrix}$where

SI means an SIMD image data set.

SP means an SIMD projection data set.

Equation (17) shows the arrangement of the SIMD data set for the uppergroup of FIG. 12. This process is referred to as packing and single SIMDinstruction data set consists of four points. This process can beapplied to the lower groups having −θ.

As noted above, the upper group and the lower group of FIG. 12 are basedon (14) and (16). Equations (14) and (16) can now be transformed intoSIMD versions as given in (17), as shown below: $\begin{matrix}{{\left( {SP} \right)_{\phi,\theta}\left( {x_{r},y_{r}} \right)} = {{\left( {SI} \right)_{x_{r}}\left( {0,y_{r}} \right)} + {\sum\limits_{n = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{\left( {1 - r_{n,\theta}} \right)\left\{ {{\left( {SI} \right)_{x_{r}}\left( {n,Z_{y_{r},n,\theta}} \right)} + {\left( {SI} \right)_{x_{r}}\left( {{- n},{Z_{y_{r},{- n},\theta} + 1}} \right)}} \right\}} +} \\{r_{n,\theta}\left\{ {{\left( {SI} \right)_{x_{r}}\left( {n,{Z_{y_{r},n,\theta} + 1}} \right)} + {\left( {SI} \right)_{x_{r}}\left( {{- n},Z_{y_{r},{- n},\theta}} \right)}} \right\}}\end{bmatrix}}}} & (18) \\{{\left( {SP} \right)_{\phi,{- \theta}}\left( {x_{r},y_{r}} \right)} = {{\left( {SI} \right)_{x_{r}}\left( {0,y_{r}} \right)} + {\sum\limits_{n = 1}^{N_{I_{x_{r}}}/2}\begin{bmatrix}{{\left( {1 - r_{n,\theta}} \right)\left\{ {{\left( {SI} \right)_{x_{r}}\left( {n,{Z_{y_{r},{- n},\theta} + 1}} \right)} + {\left( {SI} \right)_{x_{r}}\left( {{- n},Z_{y_{r},n,\theta}} \right)}} \right\}} +} \\{r_{n,\theta}\left\{ {{\left( {SI} \right)_{x_{r}}\left( {n,Z_{y_{r},{- n},\theta}} \right)} + {\left( {SI} \right)_{x_{r}}\left( {{- n},{Z_{y_{r},n,\theta} + 1}} \right)}} \right\}}\end{bmatrix}}}} & \quad\end{matrix}$

where

SI_(x) _(r) (n,Z_(y) _(r) _(,n,θ)), SI_(x) _(r) (n,Z_(y) _(r) _(,n,θ)+1)and SP_(φ,θ(x) _(r),y_(r)) are according to subgroup “a-1” of FIG. 12.

SI_(x) _(r) (−n,Z_(y) _(r) _(,−n,θ)), SI_(x) _(r) (−n,Z_(y) _(r) ,−nθ+1)and SP_(φ,θ)(x_(r),y_(r)) are according to subgroup “a-2” of FIG. 12.

SI_(x) _(r) (n,Z_(y) _(r) _(,−n,θ)), SI_(x) _(r) (n,Z_(Y) _(r)_(,−n,θ)+1) and SP_(φ,−θ)(x_(r),y_(r)) are according to subgroup “b-1”of FIG. 12.

SI_(x) _(r) (−n,Z_(r) _(r) _(,n,θ)), SI_(x) _(r) (−n,Z_(y) _(r)_(,n,θ)+1) and SP_(φ,−θ)(x_(r),y_(r)) are according to subgroup “b-2” ofFIG. 12.

Z_(y) _(r) _(,n,θ)=Z_(y) _(r) _(,−n,−θ)

Z_(y) _(r) _(m,n,θ)=Z_(y) _(r) _(,n,−θ)

This SIMD packing concept can be applied equally to a backwardprojection.

As a precaution, it should be noted that there are several special caseswithin the above symmetric pairs. These special cases indicate that somesymmetry pairs will share the same rotating image point, therebyresulting in errors. TABLE II shows these special cases, which requirespecial attention and correction.

TABLE II THE SPECIAL CASES WITHOUT SYMMETRY COUNTERPART Symmetry pairson image plane Special case φ-symmetry I(x′,y′,z_(y) _(r) _(,y′,θ)) andI(−y′,x′,z_(y) _(r) _(,y′,θ)) The points of x′ = y′ Mirror I(x′,y′,z_(y)_(r) _(,y′,θ)) and I(−x′,y′,z_(y) _(r) _(,y′,θ)) The points of x′ = 0symmetry φ-symmetry I(x′,y′,z_(y) _(r) _(,y′,θ)) and I(x′,=y′,z_(y) _(r)_(,y′,θ)) The points of y′ = 0

B. Optimization of Job Distribution in Multi-Processor Environment

The recently available PC provides a multi-processor environment. Toutilize this parallel computation capability, an instruction called“thread programming technique” can be used. The proposed SSP method iscapable of being expanded to multi-processor systems. In this case, ajob distribution with an equal amount of load is important for optimalperformance.

As the method according to the present invention is based on theassumption of circular cylindrical FOV, the job load is designed to bedistributed equally along the x′ or x_(r) axis. This is because there isuneven distribution of computational load if an equal range of x′ isassigned to each CPU, as illustrated in FIG. 6.

FIG. 6 illustrates the concept of the balanced job distribution based onthe sums of the ray path length, that is, equal distribution of thecomputational load for each group S1, S2, S3 and S4. In this case, thejob loads are divided into four groups of equal area, i.e., S1, S2, S3and S4. This redistribution suggests that each group has almost the samearea, i.e., the same amount of computational load.

Description of the Overall Method

The schematics of the new SSP projector and backprojector are providedin FIGS. 7 and 8, respectively. For SSP, the aligned frame was used byutilizing the rotation operation as well as the symmetry properties toreduce the number of loop counts. In addition, the total number ofinstructions was reduced by combining the symmetry properties and SIMDtechnique.

FIG. 7 is a flow chart showing the projection. The inner-most block 706within the loop includes z loop, the interpolation operation, theintegration along the ray path 708 and scanning along the x′ axis 710.After the projection for the entire θ range, the intermediate results ofrotated projection should be unpacked in order to be stored back to theoriginal position. It should be noted that the sizes of (I-Pack) andSIMD data sets (I-Pack+P-Pack) must be less than the cache sizes of L1and L2, respectively, in order to maximize the cache hit ratio and thedata access speed.

FIG. 8 is a flow chart showing the backprojection. First, 8 projectiondata, which have symmetric relationships, are packed for use in theSIMD. These packed projection data, P-Pack, are interpolated linearlyalong the y_(r) direction to calculate the packed backprojected imagedata, I-Pack, for the semifinal image I_(φ) (x′,y′,z) for all θdirections. As an intermediate step, these packed backprojected imagedata are unpacked and aligned to each valid position in the semifinalimage data I_(φ) (x′,y′,z) after each step of θ loop. These proceduresare repeated by loops of y_(r) and x_(r) until the semifinal image I_(φ)(x′,y′,z) is fully reconstructed. The semifinal image I_(φ) (x′,y′,z) isreconstructed for each view angle, φ. Each semi-final image I_(φ)(x′,y′,z) is rotated by its own view angle θ and then added up toreconstruct the final image I (x,y,z), the image at the originalCartesian coordinates. Similar to the projection operation, the sizes ofSIMD data sets, (I-Pack+P-Pack) and (I-Pack) should be less than thecache sizes of L2 and L1, respectively, for the maximization of thecomputation.

As shown in FIG. 7 and FIG. 8, the rotation operation was used toutilize the aligned frame. Before the projection step, the projectionplane is aligned with the aligned from (x′,y′) and each incoming imageplane is then rotated according to its respective view angle. After thebackprojection step, the original image is recovered from the alignedframe by rotating back to the original view angle, thus permitting anintermediate reconstructed image data to be formed.

To reduce processing time, the efficient use of cache was found to be animportant factor in the overall computation. Memory allocation and looporders should be optimized for the memory access pattern. Sincex_(r)(=x′) is set to the variable of the outermost loop, the memorystride on read and write can be confined to a certain range within theL2 cache size. Furthermore, the variable of the innermost loop should bez or y_(r), since most of the computation intensive operations in theSSP method are related to interpolation of z or y_(r). Therefore, it isdesirable to allocate the initial memory by the order [x′]:[y′]:[z] forbackprojection while for projection, the order should be[φ]:[x_(r)]:[θ]:[y_(r)]. These memory allocations, therefore, willdetermine the loop orders. These loop orders of projector/backprojectorare also depicted in FIG. 7 and FIG. 8. As a result of using thesesymmetries, the total number of loop counts can be reduced by half ateach step.

For the SIMD, a memory packing (step) is performed prior to eachprojection/backprojection step. After the projection step, the packedSIMD data set in the projection domain is unpacked and returned to theoriginal position. Similarly, the packed SIMD data set in the imagedomain is unpacked after each backprojection step.

Experimental Results

A. Methods

To evaluate the efficiency of the SSP method, the data or sonogram wasapplied to HRRT (High Resolution Research Tomograph developed byCPS/Siemens, Knoxville, Tenn., U.S.A). HRRT has 119,808 detectors and isdesigned for ultra high resolution brain scanning. Since HRRT has thelargest number of detectors (with the smallest detectors size) among thehuman PET scanners that are currently available, it has the highestcomputational complexity for image reconstruction. The number ofcomputations has been increased substantially (in proportion to thenumber of LORs) and is much larger compared to other existing PETscanners. The reconstruction software provided by the HRRT packagesupports iterative algorithms such as OP-OSEM3D and is the preferredreconstruction algorithm in the HRRT system. One of the advantages ofOP-OSEM3D, among others, is preventing the occurrence of negative biaswhen the random rate is high. This algorithm was chosen to test theperformance of our new SSP method since it requires the largest memorysize and is the most computation-intensive algorithm. Therefore, unlessotherwise specified, the HRRT OP-OSEM3D package S/W will be referred toas the existing method.

As can be seen, the SSP method is supported by commercially availablecomputer systems such as a common PC. Three other platforms were chosento compare the SSP method against the existing method. The firstplatform is PC 1, which has Intel dual-core 3.0 GHz, 4 GB RAM, 1 MB L2cache per CPU. The second platform, PC2, is a high-end PC with twodual-core AMD 2.4 GHz, 8 GB RAM, 1 MB L2 cache per CPU. The thirdplatform is the current HRRT computer platform, which consists of aneight node cluster system (this is denoted as HRRT CPS system), whereineach node consists of two Xeon 3.06 GHz processors with 512 KB L2 cacheand 2 GB RAM.

The concept of “span” in 3D means the mode of axial compression. Thetest was performed in span 3 and span 9, which are typically used inHRRT. The parameters for OP-OSEM3D are as follows: 256×256×207 imagepixels; six iterations; and 16 subsets. Thread programming is also used.

B. Results

Before discussing the computational speed, the accuracy of thereconstructed images (compared to the original HRRT package provided byCPS/Siemens) appears to have been emphasized. First, the results of theprojector and backprojector of the SSP were compared with the existingHRRT package. As shown in FIGS. 9 and 10, only a little difference wasfound. As mentioned previously, the SSP technique or method is basicallya tri-linear interpolation, whereas the existing method is a bi-linearmethod. FIG. 9B shows the profiles for projection data at 45° and 90°.The same comparison was made on the backprojection data as well as thefinal reconstructed images, wherein the results are shown in FIGS. 10and 11. As shown therein, the projected data and image quality of SSPare equivalent to the existing method, but yet with an advantage ofnearly two orders of magnitude in computational speed.

The performance aspect of the SSP method (i.e., execution time on PC 1and 2) was also compared. The computation time was measured for the twoexisting compression modes, namely, span 3 and span 9. The resultsindicate that the SSP method, as expected, indeed enhanced theperformance by almost 80 times, 2 (φ+90° symmetry)×2 (x_(r) ormirror-symmetry)×2 (y′-symmetry)×2 (θ-symmetry)×4 (SIMD)+otheroptimizations, compared to the existing method (original HRRT package).It should be noted that this performance enhancement can be obtainedmore or less independent of the system architecture. Furthermore, as themethod suggests, the relative improvement will be higher as more obliquerays or angles are used, as seen from TABLES III and IV.

FIG. 9A illustrates a comparison of projection data between the existingmethod and the proposed SSP method. (Top: Existing, Middle: SSP, Bottom:The difference between SSP and the existing method) FIG. 9B illustratesthe cut-views of sinogram at a specific view. (Cop: Comparison at φ=45°,Middle: Comparison at φ=90°, Bottom: The differences between SSP and theexisting method at φ=45° and φ=90°, respectively).

FIG. 10A illustrates a comparison of simple back-projection imagesbetween the existing method and the proposed SSP method. (Top: Existingmethod. Middle: New SSP technique. Bottom: The differences) FIG. 10Billustrates the cut-views (profiles) at an x axis (y=154, z=103). (Top:Comparison between the existing and new SSP. Bottom: The differencesbetween the two)

FIG. 11 illustrates a comparison of reconstructed images with sixiterations and their differences. FIG. 11A illustrates a set ofreconstructed images with the existing method (top), with new SSP method(middle) and the differences (bottom). FIG. 11B illustrates thecut-views (profiles) at an x axis (y=154, z=103).

TABLE III THE COMPARISONS OF RUNNING TIME OF THE VARIOUS OPERATIONS ONPC1. Existing new ratio Projection (span 9) 1539(s) 24(s) 64.1Backprojection (span 9) 1527(s) 25(s) 61.1 Projection (span 3) 4276(s)51(s) 83.8 Backprojection (span 3) 4255(s) 53(s) 80.3 Projection +backprojection 3066(s) 49(s) 62.6 (span 9) Projection + backprojection8531(s) 104(s)  82.0 (span 3) PC1: Intel 3.0 GHz 1 dual-core CPU, 4 GBRAM, 1 MB L2 cache per CPU.

TABLE IV THE COMPARISONS OF RUNNING TIME OF THE VARIOUS OPERATIONS ONPC2 Existing new ratio Projection (span 9)  556(s)  9(s) 61.8Backprojection (span 9)  679(s) 11(s) 61.7 Projection (span 3) 1600(s)22(s) 72.2 Backprojection (span 3) 1767(s) 22(s) 80.3 Projection +backprojection 1235(s) 20(s) 61.7 (span 9) Projection + backprojection3608(s) 44(s) 82 (span 3) PC2: AMD 2.4 Ghz 2 dual-core CPU, 8 GB RAM, 1MB L2 cache per CPU.

The use of more oblique rays or angles (e.g., span 3) results in betteraxial resolution and statistics. Since reconstruction using the EMalgorithm consists of a large number of projections and backprojections,this improvement in speed is especially important. Finally, acomparative study of the HRRT CPS system was compared with the originalHRRT CPS algorithm vs. our 2 dual core CPU (without cluster)configuration, PC2, with the SSP method and obtained a speed gain of afactor of 9 to 11. The reason for the differences between thecalculation (64 times) and the actual system based computation (9-11times) is probably due to the 16 Xeon-CPUs based original HRRT CPSsystem. The HRRT CPS system having 16 Xeon-CPUs is much more powerfulthan PC2 having only 2 dual-core CPUs. This H/W difference could havebeen the result of such reduced speed gain. In addition, SSP method caneasily extend to the cluster version.

The proposed SSP method will improve the overall computation time,especially when a dynamic functional study is desired. Currently, theimage reconstruction time for a normal HRRT PET scan with a commerciallyavailable PC (PC2 type) requires about seven to eight minutes for span 3after completing the sinogram formation and the precorrection processescompared with eighty minutes with current HRRT reconstruction package(OP-OSEM3D) utilizing an 8-node cluster system.

CONCLUSION

The present invention provides a fast method of calculating the raypath, which leads to the projection data and reconstructed image in 3Dby using the symmetry properties of the projection and image data.Exploiting these simple geometrical symmetry properties and with thehelp of SIMD, a new ultra fast reconstruction method was developed,which has a computational speed advantage of nearly two orders ofmagnitudes compared to the existing reconstruction package (specificallythat of OP-OSEM3D, the latest and most widely used PET method,especially for HRRT without compromising image quality). The keyconcepts introduced in the SSP method can be summarized as follows.First, interpolation operations are reduced from three dimension to onedimension by using rotation based projection, or the aligned frameconcept. Second, the rotation based projection is combined withsymmetric properties of the (θ,φ,y′,x_(r)) and coupled with SIMDtechnique with optimized L1 and L2 cache use. The sixteen symmetrypoints, which share the same interpolation coefficients (or theircomplementary values), were grouped, thus permitting them to beprocessed simultaneously in a projection/backprojection operation tothereby achieve a substantially improved overall reconstruction time. Tosimultaneously process these symmetric points in theprojection/backprojection, SIMD operations have also been incorporated.In particular, the SIMD scheme allowed us to simultaneously access fourdata. In addition, the data size per loop suitable for an L2 cache sizewas optimized, as well as the data structures for minimizing the memorystride.

In summary, the projection and backprojection operations were performedin the aligned frame by using the symmetry concept and SIMD with cacheoptimization and reduced the number of instructions in thereconstruction of an image from the sinogram data provided by the HRRTPET. In other words, the proposed SSP method incorporates the symmetryproperties of projection to maximum, thereby reducing the computationtime by as much as 16 times. Together with the use of the SIMD operatorand cache optimization, an overall computational speed gain by a factorclose to 80 was achieved. Such an improvement in computation time willopen the door for the dynamic functional study of high resolution PETsuch as HRRT, which suffered from an unacceptably long reconstructiontime and has been a serious undesired obstacle in molecular imaging.

This improvement will provide numerous new opportunities for PET users,especially for HRRT users in the molecular imaging community. It isclear that the proposed method helps the PET researchers to contemplatea possible dynamic study, which was limited by the unacceptable imagereconstruction time imposed by the reconstruction process. This newlydeveloped SSP method can also be applied to precorrection processes suchas attenuation and scatter correction. Finally, with the new SSP method,it will be possible to achieve true interactive scanning. The new SSPmethod can also be applied to large classes of existing PET scanners ofvarious types as well as to the cluster systems for dynamic studies.

Any reference in this specification to “one embodiment,” “anembodiment,” “example embodiment,” etc., means that a particularfeature, structure or characteristic described in connection with theembodiment is included in at least one embodiment of the invention. Theappearances of such phrases in various places in the specification arenot necessarily all referring to the same embodiment. Further, when aparticular feature, structure or characteristic is described inconnection with any embodiment, it is submitted that it is within thepurview of one skilled in the art to effect such feature, structure orcharacteristic in connection with other ones of the embodiments.

Although embodiments have been described with reference to a number ofillustrative embodiments thereof, it should be understood that numerousother modifications and embodiments can be devised by those skilled inthe art that will fall within the spirit and scope of the principles ofthis disclosure. More particularly, various variations and modificationsare possible in the component parts and/or arrangements of the subjectcombination arrangement within the scope of the disclosure, the drawingsand the appended claims. In addition to variations and modifications inthe component parts and/or arrangements, alternative uses will also beapparent to those skilled in the art.

1. A method for reconstructing 3D image, comprising: detecting aplurality of line of responses (LORs) emitted from an object;transforming said plurality of LORs into first sinogram data;back-projecting said first sinogram data with a plurality of projectionangles to produce image data for said object; and projecting saidproduced image data with said plurality of projection angles totransform said image data into second sinogram data, wherein saidback-projecting includes filling pixels of image plane for each of saidplurality of projection angles with said first sinogram data androtating a coordinate axis of said image plane with correspondingprojection angle to produce said image data, wherein said projectingincludes rotating said image data with corresponding projection angle inan opposite direction before projecting said image data with saidplurality of projection angles, and wherein said projecting and saidback-projecting use symmetry properties in coordinate space.
 2. Themethod of claim 1, wherein said symmetry properties include at least oneof symmetry properties in view angles (φ-symmetry), symmetry propertiesin oblique angles (θ-symmetry) and symmetry properties in interpolationcoefficients along the projection lines (y′-symmetry).
 3. The method ofclaim 2, wherein said projecting and back-projecting include executing aplurality of projecting operations having said symmetry propertiessimultaneously by using single instruction multiple data (SIMD).
 4. Themethod of claim 3, wherein said using SIMD includes arranging data sethaving same oblique angle in amplitude.
 5. The method of claim 2,wherein said projecting and back-projecting include dividing said imagedata into a plurality of equal subsets, and at each of a plurality ofprocessors, performing projecting operation for each of said pluralityof subsets.
 6. The method of claim 5, wherein said dividing said imagedata into a plurality of equal subsets includes distributing jobs sothat sums of ray path length is equal.
 7. The method of claim 1, whereinsaid projecting and said back-projecting include executing a pluralityof projecting operation having said symmetry properties simultaneouslyby using single instruction multiple data (SIMD).
 8. The method of claim7, wherein said using SIMD includes arranging data set having sameoblique angle in amplitude.
 9. The method of claim 1, wherein saidprojecting and back-projecting include dividing said image data into aplurality of equal subsets, and at each of a plurality of processors,performing projecting operation for each of said plurality of subsets.10. The method of claim 9, wherein said dividing said image data into aplurality of equal subsets includes distributing jobs so that sums ofray path length is equal.
 11. An apparatus for reconstructing 3D image,comprising: means for detecting a plurality of line of responses (LORs)emitted from an object; means for transforming said plurality of LORsinto first sinogram data; means for back-projecting said first sinogramdata with a plurality of projection angles to produce image data forsaid object; and means for projecting said produced image data with saidplurality of projection angles to transform said image data into secondsinogram data, wherein said back-projecting includes filling pixels ofimage plane for each of said plurality of projection angles with saidfirst sinogram data and rotating a coordinate axis of said image planewith corresponding projection angle to produce said image data, whereinsaid projecting includes rotating said image data with correspondingprojection angle in an opposite direction before projecting said imagedata with said plurality of projection angles, and wherein saidprojecting and said back-projecting use symmetry properties incoordinate space.
 12. The apparatus of claim 11, wherein said symmetryproperties include at least one of symmetry properties in view angles(φ-symmetry), symmetry properties in oblique angles (θ-symmetry) andsymmetry properties in interpolation coefficients along the projectionlines (y′-symmetry).
 13. The apparatus of claim 12, wherein saidprojecting and said back-projecting include executing a plurality ofprojecting operation having said symmetry properties simultaneously byusing single instruction multiple data (SIMD).
 14. The apparatus ofclaim 13, wherein said using SIMD includes arranging data set havingsame oblique angle in amplitude.
 15. The apparatus of claim 12, whereinsaid projecting and said back-projecting include dividing said imagedata into a plurality of equal subsets, and at each of a plurality ofprocessors, performing projecting operation for each of said pluralityof subsets.
 16. The apparatus of claim 15, wherein said dividing saidimage data into a plurality of equal subsets includes distributing jobsso that sums of ray path length is equal.
 17. The apparatus of claim 11,wherein said projecting and back-projecting include executing aplurality of projecting operations having said symmetry propertiessimultaneously by using single instruction multiple data (SIMD).
 18. Theapparatus of claim 17, wherein said using SIMD includes arranging dataset having same oblique angle in amplitude.
 19. The apparatus of claim11, wherein said projecting and back-projecting include dividing saidimage data into a plurality of equal subsets, and at each of a pluralityof processors, performing projecting operation for each of saidplurality of subsets.
 20. The apparatus of claim 19, wherein saiddividing said image data into a plurality of equal subsets includesdistributing jobs so that sums of ray path length is equal.
 21. Anapparatus for reconstructing a 3D image, the apparatus comprising:apparatus for detecting a plurality of line of responses (LORs) emittedfrom an object; and a processor programmed with an algorithm for:transforming the plurality of LORs into first sinogram data;back-projecting the first sinogram data with a plurality of projectionangles to produce image data for the object, the back-projectingincluding filling pixels of an image plane for each of said plurality ofprojection angles with said first sinogram data and rotating acoordinate axis of the image plane with a corresponding projection angleto produce the image data; and projecting the image data with theplurality of projection angles to transform the image data into secondsinogram data, the projecting including rotating the image data with acorresponding projection angle in an opposite direction beforeprojecting the image data with the plurality of projection angles,wherein the projecting and back-projecting use symmetry properties incoordinate space.